For more details on using the NetCoMi package see https://github.com/stefpeschel/NetCoMi.
## Warning: package 'phyloseq' was built under R version 3.6.1
## [1] '1.30.0'
The following code was run on a cluster.
# # load phyloseq data
# load('/n/home05/ajsommer/NetCoMi_cluster/ps_to_net.RData')
# gut_split <- metagMisc::phyloseq_sep_variable(ps_prune, "W")
# net_W <- netConstruct(gut_split$`0`, gut_split$`1`, verbose = 2,
# filtTax = "highestVar",
# filtTaxPar = list(highestVar = 50),
# measure = "spieceasi",
# measurePar = list(method = "glasso",
# nlambda=20,
# pulsar.params=list(rep.num=50)),
# normMethod = "none", zeroMethod = "none",
# sparsMethod = "none", seed = 123456, matchDesign = c(1,1))
# props_W <- netAnalyze(net_W, clustMethod = "cluster_fast_greedy")
We plot the results of the netAnalyze function.
phyl_ps_prune <- as.factor(tax_table(ps_prune)[, 'Phylum'])
names(phyl_ps_prune) <- rownames(tax_table(ps_prune)[, 'Phylum'])
phylcol <- rainbow(length(unique(phyl_ps_prune)))
set.seed(1)
props_W_igraph <- graph_from_adjacency_matrix(props_W$input$assoMat1,
#mode = "undirected",
weighted = TRUE)
FA_coord <- layout.forceatlas2(props_W_igraph, iterations=2000, plotstep=2001)
rownames(FA_coord) <- names(phyl_ps_prune)
plot(props_W,
layout = FA_coord,
# layout = "layout_with_fr",
sameLayout = TRUE, layoutGroup = 1, labels = FALSE,
featVecCol = phyl_ps_prune, nodeColor = "feature",nodeTransp = 55,
borderCol = "gray40", highlightHubs = FALSE,nodeSize = "normCounts",
nodeSizeSpread = 1, cexNodes = 5, cexTitle = .8,
edgeTranspLow = 70, edgeTranspHigh = 50,
mar = c(3,5,3,2), groupNames = c("Polluted air", "Clean air"))
legend("center", inset = c(-1,-10),
legend=levels(phyl_ps_prune),
col=phylcol, bty="n", text.col = phylcol, cex = .5)
legend("bottom", inset = c(0,-2),
legend=c("positive","negative"), lty = 1,
col=c("darkgreen","red"), bty="n", cex =.6)
nclust <- max(as.numeric(names(table(props_W$clustering$clust1))))+1
col <- polychrome(nclust)
names(col) <- 1:nclust
clusters <- as.factor(c(props_W$clustering$clust1,props_W$clustering$clust2))
col_assign <- col[clusters]
names(col_assign) <- names(clusters)
# which(tax_table(ps_prune)[,"Genus"] == "Blautia")
# tax_table(ps_prune)[which(tax_table(ps_prune)[,"Genus"] == "Bacteroides"),]
# 248
shapeVec <- rep(1, ncol(otu_table(ps_prune)))
shapeVec[248] <- 2
names(shapeVec) <- names(phyl_ps_prune)
plot(props_W,
layout = FA_coord,
# layout = "layout_with_fr",
sameLayout = TRUE, layoutGroup = 1,
labels = list(props_W$clustering$clust1, props_W$clustering$clust2),
labelFont = 1, cexLabels = 2,
nodeColor = "cluster",
colorVec = col,
nodeTransp = 40,
borderCol = "gray40", highlightHubs = FALSE,
nodeSize = "fix", cexNodes = 3,
nodeShape = c("circle", "triangle"),
featVecShape = shapeVec,
edgeTranspLow = 80, edgeTranspHigh = 50,
groupNames = c("Polluted air", "Clean air"))
names_clust <- as.character(names(which(props_W$clustering$clust1 == 7)))
labels_phyl <- substr(as.character(phyl_ps_prune[names_clust]), 1, 3)
plot(props_W,
# layout = FA_coord,
sameLayout = TRUE, layoutGroup = 1,
labels = labels_phyl,
labelFont = 1, cexLabels = 1.5,
featVecCol = phyl_ps_prune[names_clust],
nodeColor = "cluster",
colorVec = col,
nodeTransp = 40,
borderCol = "gray40", highlightHubs = FALSE,
nodeSize = "fix", cexNodes = 1,
nodeFilter = "names",
nodeShape = c("circle", "triangle"),
featVecShape = shapeVec,
nodeFilterPar = names_clust,
edgeTranspLow = 80, edgeTranspHigh = 50,
groupNames = c("Polluted air", "Clean air"))
col1 <- colorRampPalette(c("#7F0000", "red", "#FF7F00", "yellow", "white",
"cyan", "#007FFF", "blue", "#00007F"))
cor_0 <- props_W$input$assoMat1[names_clust,names_clust]
cor_1 <- props_W$input$assoMat2[names_clust,names_clust]
par(mfrow = c(1,2))
corrplot(cor_0, method = "color", tl.col = "grey", mar=c(0,0,5,0),
col = col1(50), title = "Polluted air")
corrplot(cor_1, method = "color", tl.col = "grey", mar=c(0,0,5,0),
col = col1(100), title = "Clean air")